clear all
set more off
set matsize 10000
set maxvar 10000
adopath + ../code/gslab_tools/

preliminaries, doutf(../figures)
graph set window fontface "Times New Roman"

****************
** Figure 1  ****
****************

*Figure 1 does not use data


****************
** Figure 2 ****
****************

use "../derived/Combined/combined_arm_level.dta" if year!=., clear
gen relative_year = year - approve_year
drop if relative_year==.
drop if relative_year>44                                                 
gen n=1
collapse (mean) sponsor (sum) n, by(relative_year)

twoway (scatter sponsor relative_year, yaxis(1) mcolor(gs6) msymbol(diamond)), ///
	ytitle("Share of Arms Sponsored", axis(1)) xtitle("Years since FDA Approval") ////
	scheme(s2mono) graphregion(color(white))
graph export ../figures/Figure2a.eps, replace

twoway (bar n relative_year, yaxis(2) fintensity(10) ), ///
	ytitle("# of Arms", axis(2)) xtitle("Years since FDA Approval") ///
	scheme(s2mono) graphregion(color(white))
graph export ../figures/Figure2b.eps, replace



****************
** Figure 3 ****
****************

*Antidepressant Panel
use ../derived/Diff_in_diff_setup/two_by_two_w_mean.dta, clear
replace drug="All Active vs. Placebo" if drug=="" 
gen first= drug=="All Active vs. Placebo"
sort first dd

gen label_drugs=_n
count
forvalues x = 1/`r(N)'{
label define label_drugs `x' "`=drug[`x']'", add
}
label values label_drugs label_drugs

graph twoway (scatter label_drugs share_respond3, mc(black)) ///
	(scatter label_drugs share_respond4, mc(black) ms(Oh)) ///
	(scatter label_drugs share_respond1, msymbol(triangle) mc(gs9)) ///
	(scatter label_drugs share_respond2, mc(gs9) ms(th)), ///
	ylabel(1(1)6, nolab) xtitle("Share Respond") ytitle("") ///
	legend(label(1 "Drug, Sponsored") label(2 "Placebo, Sponsored") ///
	label(3 "Drug, Not Sponsored") label(4 "Placebo, Not Sponsored")  ///
	region(lstyle(none) col(white)) nobox) scheme(s1mono) ///
	plotregion(color(white)) graphregion(color(white) lcolor(white)) ///
	bgcolor(white) saving(points, replace)

graph twoway (scatter label_drugs dd, mc(white)) ///
	(scatter label_drugs dd, mc(gs2)), xline(0, lpattern(dash) lcolor(gs10)) ///
	ylabel(1(1)6, angle(horizontal) valuelabel) ///
	xtitle("Sponsorship Effect", margin(0 0 0 2.7)) ytitle("") ///
	legend(order(2 1) col(1) label(2 "Diff-in-Diff") label(1 "") ///
	symysize(20) region(lstyle(none) col(white) lcolor(white)) nobox) ///
	plotregion(color(white)) graphregion(color(white) lcolor(white)) ///
	bgcolor(white) saving(dd, replace)

graph combine dd.gph points.gph, plotregion(color(white)) graphregion(color(white) lcolor(white) )
graph export ../figures/Figure3a.eps, replace


use ../derived/Diff_in_diff_setup/two_by_two_active_anti_w_mean.dta, clear
replace drug="All Active vs. Active" if drug==""
gen first= drug=="All Active vs. Active"
sort first dd

gen label_drugs=_n
count
forvalues x = 1/`r(N)'{
label define label_drugs `x' "`=drug[`x']'", add
}
label values label_drugs label_drugs

graph twoway (scatter label_drugs dd, mc(white)) ///
	(scatter label_drugs dd, mc(gs2)), title("Antidepressants") ///
	xline(0, lpattern(dash) lcolor(gs10)) ///
	ylabel(1(1)26, angle(horizontal) valuelabel) ///
	xtitle("Sponsorship Effect") ytitle("") ///
	legend(order(2 1) col(1) label(2 "Diff-in-Diff") label(1 "") ///
	symysize(20) region(lstyle(none) col(white) lcolor(white)) nobox) ///
	plotregion(color(white)) graphregion(color(white) lcolor(white)) ////
	bgcolor(white) saving(dd_anti, replace)

*Antipsychotic Panel
use ../derived/Diff_in_diff_setup/two_by_two_active_schiz_w_mean.dta, clear
replace drug="All Active vs. Active" if drug==""
gen first= drug=="All Active vs. Active"
sort first dd

gen label_drugs=_n
count
forvalues x = 1/`r(N)'{
label define label_drugs `x' "`=drug[`x']'", add
}
label values label_drugs label_drugs

graph twoway (scatter label_drugs dd, mc(white)) ///
	(scatter label_drugs dd, mc(gs2)), title("Antipsychotics") ///
	xline(0, lpattern(dash) lcolor(gs10)) ///
	ylabel(1(1)13, angle(horizontal) valuelabel) ///
	xtitle("Sponsorship Effect") ytitle("") ///
	legend(order(2 1) col(1) label(2 "Diff-in-Diff") label(1 "") ////
	symysize(20) region(lstyle(none) col(white) lcolor(white)) nobox) ///
	plotregion(color(white)) graphregion(color(white) lcolor(white)) ///
	bgcolor(white) saving(dd_schiz, replace)

graph combine dd_anti.gph dd_schiz.gph, plotregion(color(white)) graphregion(color(white) lcolor(white) )
graph export ../figures/Figure3b.eps, replace

****************
** Figure 4 ****
****************

use "../derived/Combined/combined_pair_level.dta", clear
gen published = year!=.

reg published stdz_outcome_relative_all [aw=weight] if sponsor ==1
local sponsor_beta = round(_b[stdz_outcome_relative_all], 0.001)
local sponsor_se = round(_se[stdz_outcome_relative_all], 0.001)
reg published stdz_outcome_relative_all [aw=weight] if sponsor ==0
local nosponsor_beta = round(_b[stdz_outcome_relative_all], 0.001)
local nosponsor_se = round(_se[stdz_outcome_relative_all], 0.001)

xtile efficacy_bin=stdz_outcome_relative_all [aw=weight], nq(25)


collapse (mean) published stdz_outcome_relative_all (sum) weight [aw=weight], by(efficacy_bin sponsor)
twoway (scatter published stdz_outcome_relative_all if sponsor==1, mcolor(black)) ///
	(scatter published stdz_outcome_relative_all if sponsor==0, msymbol(plus) mcolor(gs9)) ///
	(lfit published stdz_outcome_relative_all if sponsor==1, lcolor(black) lpattern(dash)) ///
	(lfit published stdz_outcome_relative_all if sponsor==0, lcolor(gs9) lpattern(shortdash_dot)), ////
	legend(order(1 "Sponsored" 2 "Not Sponsored" 3 "`sponsor_beta' (`sponsor_se')" ///
	4 "`nosponsor_beta' (`nosponsor_se')")) scheme(s2mono) ///
	xtitle("Relative Efficacy (Standardized)") ytitle("Share Published") ///
	graphregion(color(white)) 
graph export ../figures/Figure4.eps, replace


****************
** Figure 5 ****
****************

use "../derived/NCT_registry/registry_link.dta", clear
keep if year!=.
gen n=1
drop if year<1980
collapse (sum) n (mean) linked, by(year)

twoway 	(bar n year, yaxis(2) fintensity(10) ytitle("# of Trial Arms", axis(2))) ///
	(scatter linked year, yaxis(1) ylabel(0(0.2)1) ///
	mcolor(gs6) ytitle("Share Pre-Registered", axis(1))) ////
	(pci 0 2005.4 32 2005.4, yaxis(2) lpattern(dash) lcolor(gs6)), ///
	legend(order(1 "# of Trial Arms" 3 "Share Pre-Registered")) ///
	scheme(s2mono) graphregion(color(white)) xlabel(1980(5)2015)
graph export ../figures/Figure5a.eps, replace
	
	
use "../derived/Combined/combined_pair_level.dta", clear
merge m:1 studyname using "../derived/NCT_registry/registry_link", ///
	keepusing(linked) nogen
egen drug_group=group(drug)
drop if year==.
		
drop year_group
gen year_group1 = year<=1991
gen year_group2 = year >=1992 & year<=1993
gen year_group3 = year >=1994 & year<=1995
gen year_group4 = year >=1996 & year<=1997
gen year_group5 = year >=1998 & year<=1999
gen year_group6 = year >=2000 & year<=2001
gen year_group7 = year >=2002 & year<=2003
gen year_group8 = year >=2004 & year<=2005
gen year_group9 = year >=2006 & year<=2007
gen year_group10 = year >=2008 & year<=2009
gen year_group11 = year >=2010 & year<=2011
gen year_group12 = year >=2012 & year<=2013
gen year_group13 = year >=2014 & year!=.
gen year_group14 = year==.
gen post = year>=2006 & year!=.

forvalues yr = 1/13{
	gen sponsor_`yr' = sponsor * year_group`yr'
}
gen sponsormi = sponsor*year_group14
egen sponsor_0 = rowtotal(sponsor_1 - sponsor_8)

label var sponsor_1 "1991 & Before"
label var sponsor_2 "1992-1993"
label var sponsor_3 "1994-1995"
label var sponsor_4 "1996-1997"
label var sponsor_5 "1998-1999"
label var sponsor_6 "2000-2001"
label var sponsor_7 "2002-2003"
label var sponsor_8 "2004-2005"
label var sponsor_9 "2006-2007"
label var sponsor_10 "2008-2009"
label var sponsor_11 "2010-2011"
label var sponsor_12 "2012-2013"
label var sponsor_13 "2014 & After"

areg new_relative sponsor_1-sponsor_13 sponsormi o.sponsor i.scale_group ///
	year_group* [aw=weight], a(factor) cluster(cluster_group)
estimates store full

coefplot full, keep(sponsor_*) vertical omitted scheme(s2mono)  ///
	xtitle("Publication Year") ytitle("Sponsorship Effect") ///
	coeflabels(,angle(45)) yline(0) xline(8.5, lpattern(dash) lcolor(gs8)) ///
	graphregion(color(white)) ylabel(,nogrid)
graph export ../figures/Figure5b.eps, replace	


****************
** Figure 6 ****
****************

use "../derived/Combined/combined_pair_level.dta" if year!=., clear
merge m:1 drug using "../derived/MEPS/meps_spending", assert(1 3) ///
	nogen keepusing(first_five_n)
egen drug_group=group(drug)
areg new_relative o1.sponsor##i.drug_group i.scale_group i.year_group [aw=weight], ///
		a(factor) cluster(cluster_group)

gen coeff = .
gen base_coeff=.
levelsof drug_group, local(drug_levels)
foreach i of local drug_levels {
	replace coeff = _b[1.sponsor#`i'.drug_group]+_b[1.sponsor] if drug_group ==`i'
	replace base_coeff = _b[1.sponsor] if drug_group ==`i'
}

gen base_stdz_outcome = stdz_outcome_all if sponsor==0

gen n=1
collapse (sum) weight (first) coeff base_coeff approve* firm* drug_type* ///
	(mean) base_stdz_outcome, by(drug drug_group *five* patent*)	
	
sort coeff
gen pos =3
save ../derived/Combined/combined_pair_level_meps_merge.dta, replace

bys drug_type (approve_year): gen order=_n
gen first_inclass = order==1

gen class_year_temp = approve_year if first_inclass
bys drug_type: egen class_year=min(class_year_temp)
drop *temp
gen rel_class_year = approve_year - class_year

drop if coeff==base_coeff 
sort approve_year
gen overall_order=_n

drop if approve_year==.

replace pos = 3 if drug=="citalopram"
replace pos = 12 if drug=="fluvoxamine"
replace pos = 6 if drug=="fluoxetine"
replace pos = 3 if drug=="mirtazapine"
replace pos = 5 if drug=="sertraline"
replace pos = 1 if drug=="amitriptyline"
replace pos = 12 if drug=="risperidone"
replace pos = 9 if drug=="bupropion"
replace pos = 6 if drug=="duloxetine"
replace pos = 6 if drug=="trazodone"
replace pos = 4 if drug=="quetiapine"
replace pos = 4 if drug=="venlafaxine"
replace pos = 9 if drug=="escitalopram"
replace pos = 6 if drug=="venlafaxine"
replace pos = 1 if drug=="ziprasidone"
replace pos = 1 if drug=="aripiprazole"
replace pos=6 if drug=="clozapine"
replace pos=6 if drug=="olanzapine"
replace pos=10 if drug=="paroxetine"

reg coeff rel_class_year [aw=weight]
local beta = round(_b[rel_class_year], 0.001)
local se = round(_se[rel_class_year], 0.001)

twoway	(scatter coeff rel_class_year, mlabel(drug) mlabv(pos) msymbol(none) mlabg(3)) ///
	(scatter coeff rel_class_year [w=weight], msymbol(circle_hollow)) ///
	(lfit coeff rel_class_year [aw=weight], lpattern(dash)), ytitle("Sponsorship Coefficient") ///
	scheme(s1mono) xtitle("Year Since First Drug Approved in Class") xscale(range(-1 15)) ///
	legend(off) note("Best Fit Line Slope: `beta' (`se'), P=0.026")
graph export ../figures/Figure6.eps, replace


****************
** Figure 7 ****
****************

use ../derived/Combined/combined_pair_level_meps_merge.dta, clear

replace pos = 3 if drug=="citalopram"
replace pos = 12 if drug=="fluvoxamine"
replace pos = 9 if drug=="fluoxetine"
replace pos = 3 if drug=="mirtazapine"
replace pos = 5 if drug=="sertraline"
replace pos = 9 if drug=="amitriptyline"
replace pos = 11 if drug=="risperidone"
replace pos = 3 if drug=="bupropion"
replace pos = 6 if drug=="duloxetine"
replace pos = 6 if drug=="trazodone"
replace pos = 4 if drug=="quetiapine"
replace pos = 4 if drug=="venlafaxine"
replace pos = 9 if drug=="escitalopram"
replace pos = 2 if drug=="venlafaxine"
replace pos = 6 if drug=="ziprasidone"
replace pos = 4 if drug=="olanzapine"

drop if coeff==base_coeff // Drop if sponsorship effect not estimated within this drug 

reg coeff first_five_n [aw=weight]
local beta = round(_b[first_five_n], 0.000001)
local se = round(_se[first_five_n], 0.000001)

twoway (scatter coeff first_five_n, mlabel(drug) mlabv(pos) msymbol(none) mlabg(3)) ///
	(scatter coeff first_five_n [w=weight], msymbol(circle_hollow)) ///
	(lfit coeff first_five_n [aw=weight], lpattern(dash)), ytitle("Sponsorship Coefficient") ///
	scheme(s1mono) xtitle("Average MEPS Prescriptions Post Approval") ///
	legend(off) note("Best Fit Line Slope: `beta' (`se'), P=0.057")
graph export ../figures/Figure7.eps, replace


cap erase points.gph
cap erase dd.gph
cap erase dd_anti.gph
cap erase dd_schiz.gph
